home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 188_01 / trans.doc < prev    next >
Text File  |  1987-09-29  |  14KB  |  275 lines
  1. .mt 0
  2. .mb 11
  3. .po 13
  4. /*
  5. HEADER:        ;
  6. TITLE:        C elementary transcendentals;
  7. VERSION:    1.0;
  8.  
  9. DESCRIPTION:       Source    code    for   all   standard    C                            
  10.                transcendentals
  11.  
  12.         Employs  ldexp()  and  frexp()  functions;  if                
  13.                suitable  versions  of these are  not  provided                
  14.                by  a given compiler,  the versions provided in                
  15.                source  code  will require  adaptation  to  the                
  16.                double float formats of the compiler.
  17.  
  18. KEYWORDS:    Transcendentals, library;
  19. SYSTEM:        CP/M-80, V3.1;
  20. FILENAME:    TRANS.C;
  21. WARNINGS:      frexp()   and  ldexp()  are   implementation                
  22.                dependent.   The  compiler  employed  does  not 
  23.                support    minus   (-)   unary   operators   in                
  24.                initializer  lists,  which are required by  the                
  25.                code.
  26. AUTHORS:    Tim Prince;
  27. COMPILERS:    MIX C, v2.0.1;
  28. */
  29.               Transcendental Function Algorithms
  30.  
  31.      All  programming languages which are descended from  FOR-
  32. TRAN  or  Algol include pre-defined transcendental  functions.     
  33. Standard  textbooks on programming assume that the usual  math 
  34. functions  are already available and not worth the  paper  re-
  35. quired  to describe how they might be coded.   This is  unfor-
  36. tunate,  as few computer systems provide reliable transcenden-
  37. tals,  and  a basic understanding of their contents is helpful 
  38. in their use.
  39.  
  40.      The  following  examples use polynomial  representations, 
  41. because  they lead to compact code and give  average  accuracy 
  42. within one bit of the best available.  Faster algorithms exist 
  43. but require more code or data storage,  if written in the same 
  44. language.   Commercial  versions of the functions are  usually 
  45. written  in  assembler,  but the ideas are best  expressed  in 
  46. higher  level languages.   Obscure tricks which  an  anonymous 
  47. coder thought would improve speed are useless when the program 
  48. fails  and the source code is not available.   These functions 
  49. have been tested in nearly the same form on several  computers 
  50. in both FORTRAN and C.
  51.  
  52.      The  author's hope is that these examples may set a mini-
  53. mum standard for accuracy and speed of C library  transcenden-
  54. tals.   Failing  this,  the  individual user should  have  the 
  55. opportunity  to check whether any deficiencies of his code may 
  56. be due to the library provided with the C compiler.  C's  soon 
  57. to  be standardized support of bit fields and low level opera-
  58. tions on doubles helps in the expression of these  algorithms.  
  59. Since  C,  even in the future ANSI standard,  will not support 
  60. pure single precision on computers which were not designed for 
  61. mixed precision arithmetic, and the standard library functions 
  62. are double, we show only the double precision, with the under-
  63. standing that changes in the polynomials would adapt to  other 
  64. precisions.   The file TRANSLIB.FOR shows FORTRAN versions  in 
  65. single precision, which should serve as a demonstration of the 
  66. changes required to vary precision and language.
  67.  
  68.      The  polynomial  coefficients were determined by  running 
  69. one  of  Steve Ruzinsky's programs (1).   In  each  case,  the 
  70. polynomial  is derived by truncating a standard Taylor  series 
  71. or  rational  polynomial representation of  the  function  and 
  72. using Steve's program to obtain more accuracy with fewer terms 
  73. in the required range.   In some cases, 25 digit arithmetic is 
  74. required  to fit a polynomial accurate to at least 16  digits.  
  75. Tables of coefficients have been generated for polynomials  of 
  76. any accuracy required up to 20 or more significant digits.
  77.  
  78.      Trigonometric  functions can be calculated with  adequate 
  79. efficiency  using portable code,  even in archaic dialects  of 
  80. FORTRAN,  Pascal,  and C.   When using polynomials, it appears 
  81. best  to use the sin() function as the basic series and calcu-
  82. late the others from it.   A polynomial for tan() requires  as 
  83. many terms as sin() and cos() together.  Since cos() and tan() 
  84. may be calculated from sin(),  only one shorter table of coef-
  85. ficients  is  needed.   In effect,  tan() is calculated  by  a 
  86. rational polynomial approximation.   Expressing these trigono-
  87. metric relationships by #define teaches the preprocessor alge-
  88. bra  which  an  optimizing compiler could  use  to  advantage.  
  89. Check  that your compiler does not execute functions  multiple 
  90. times  as a result of macro expansion.   Tan() could be calcu-
  91. lated much faster with a big table as is built in to the 8087, 
  92. but  most programs we have seen spend more time on  sin()  and 
  93. cos().
  94.  
  95.      The  first  step in calculating sin() is  to  reduce  the 
  96. range to |x| < pi/2 by subtracting the nearest multiple of pi.  
  97. The  ODD() function shown assumes that (int) will extract  the 
  98. low  order  bits  in case the argument  exceeds  maxint.   The 
  99. multiplication  and subtraction should be done in long  double 
  100. if it is available, but in any case the range reduction should 
  101. not reduce accuracy when the original argument was already  in 
  102. the  desired  range.   The  method most often used  saves  one 
  103. divide (which may be changed to a multiply), at the cost of an 
  104. unnecessary  roundoff.   Worse,  the result may  underflow  to 
  105. zero.   Underflows  in the code shown occur only where a  zero 
  106. does not affect the result.  With the precautions taken, extra 
  107. precision   is  not  needed  anywhere  other  than  the  range 
  108. reduction.
  109.  
  110.      Most  computer  systems include a  polynomial  evaluation 
  111. function,  sometimes in microcode.  To save space we prefer to 
  112. use such a function even though it may be slower than  writing 
  113. out  the  polynomial in Horner form.   Such functions  operate 
  114. like poly() shown here,  but (presumably) faster.   This modu-
  115. larity  permits  special  hardware features to  be  used  more 
  116. effectively,  although the time required to call another func-
  117. tion may be excessive.  Loop unrolling (shown here) or odd and 
  118. even  summing (see exp2()) are often  effective.   The  series 
  119. used  for elementary functions have terms increasing in magni-
  120. tude fast enough so that there is no need for extra  precision 
  121. intermediate computation.
  122.  
  123.      The  arctan function uses several substitutions to reduce 
  124. to a range of 0 < x < 1.   Although the Taylor series does not 
  125. converge over this range,  the minimax polynomial does, but so 
  126. many terms are needed that one further step of range reduction 
  127. is preferred.   Code length may be traded for speed by using a 
  128. table of tangents with a shorter series.   A rational  polyno-
  129. mial(3)  is used by many run-time libraries,  at some cost  of 
  130. accuracy.   A  series with one less term would still give over 
  131. 15  digits  accuracy,  which is the maximum possible  on  some 
  132. systems.
  133.  
  134.      Although  log base 2 is not usually included in the  list 
  135. of  library functions,  it is almost always used as the  basic 
  136. logarithm  for binary floating point,  because it is the  most 
  137. efficient  base for exponentiation with a float  power.   Many 
  138. implementations fail to obtain the logarithm of numbers  close 
  139. to 1.   The following routine is the most widely tested of the 
  140. examples  given here,  and has been found accurate on all  ma-
  141. chines which subtract 1 accurately (some don't!).   In theory, 
  142. precision  is  lost when the exponent is added in at the  end, 
  143. but  this  is not serious in the most common cases of  numbers 
  144. near 1.  Making the function long double would correct this.
  145.  
  146.      Log2  can  be written in portable code  if  the  standard 
  147. functions  frexp() and ldexp() for extracting and changing the 
  148. exponent field of a floating point number are  available.   On 
  149. many  systems,  adequate speed will not be obtained unless the 
  150. code  is  supplied in line.   This should  be  possible  using 
  151. unions of structures,  but direct use of any applicable float-
  152. ing point hardware ins